## This file is for analysing spatial panel data set "tx_emp_wind.csv" using R package "splm"
## Clean the R memory
rm(list = ls())

## Load library packages
library("spdep");         # Define listw class, and transform it into a sparse matrix
library("plm")            # Panel data model
library("splm");          # Spatial panel data model
library("xtable");        # Export results to Latex
library("foreign");       # Inport data from Stata,etc.
library("calibrate")      # Show labels to specific points in plots

## Export the following output to a txt file.
sink(file = "~/HOLA/shale_job/txWind.txt", append = FALSE, type = c("output", "message"), split = FALSE)

## Import data as a data.frame.
txData <- read.csv("~/HOLA/research/green-jobs/data/tx_emp_windcumu.csv", header = T) 
## load("~/HOLA/shale_job/txData.RData")
txData$RRC = factor(txData$RRC);
txData$windInd = factor(txData$windInd);
txData$windInd[txData$cumuCap > 0] = 1;
## Import spatial weights matrix.
wmat = as.matrix(read.csv("~/HOLA/research/green-jobs/data/wmat/adjmatrix.csv", header = TRUE, row.names = 1))

## Look over the data set
names(txData);
dim(txData);
head(txData);
tail(txData);
summary(txData);

## Create panel data frame
txPanel <- pdata.frame(txData, index = c("county", "year"), drop.index = TRUE, row.names = TRUE)

pdf("~/HOLA/research/green-jobs/latex/figures/tx_wind_pair.pdf", height=6, width=8)
pairs(txPanel)
dev.off()

## summary indicates the total variation of the variable and the share of this variation that is due to the individual and the time dimensions.
summary(txPanel$emp)
summary(txPanel$capacity)
## as.matrix indicates the total variation of the variable and the share of this variation that is due to the individual and the time dimensions. 
head(as.matrix(txPanel$capacity));
## head(lag(txPanel$emp, 0:2))
## head(diff(txPanel$emp), 15)
head(Within(txPanel$emp),15)
head(between(txPanel$emp),15)
## Basic regression model
fm <- emp ~ + windInd - capacity - cumuCap

## Pooling effects:

job.pool <- plm(fm,txPanel,model = "pooling")
summary(job.pool)

## Fixed effects:

## One-way
job.fe <- plm(fm,txPanel,model = "within")
summary(job.fe)

## Two-way
job.twfe <- plm(fm,txPanel,model = "within",effect = "twoways")
summary(job.twfe);
fe.spatial = fixef(job.twfe, type = "dmean");

## pdf("~/HOLA/shale_job/figures/spatialEf.pdf", height=6, width=8)
plot(fe.spatial,main = "County-specific fixed effects", xlab = "county", ylab = "Fixed effects")
textxy(fe.spatial, names(fe.spatial))
## dev.off()

fe.time = fixef(job.twfe, effect = "time");

## pdf("~/HOLA/shale_job/figures/timeEf.pdf", height=6, width=8)
plot(fe.time,main = "Time effects", xlab = "Year", ylab = "Time effects", type = 'o', lwd = 1)
## dev.off()

## Random effects:

## One-way effect, from Swamy and Arora (1972)
job.re <- plm(fm,txPanel,model = "random")
summary(job.re)
## Two-way effects, from Amemiya (1971),
job.twre <- plm(fm,txPanel,model = "random",random.method = "amemiya", effect = "twoways")
summary(job.twre)
## attach(job.twre)
## pdf("~/HOLA/research/green-jobs/latex/figures/tx_wind_twre.pdf", height=6, width=6)
## plot(txPanel$wells,txPanel$emp, main = "Random effects model")
## abline(coefficients, lwd = 2, col = "blue")
## mtext(bquote( emp == .(coefficients[2]) * txPanel$wells + .(coefficients[1])), side=3, line=0) 
## dev.off()
## detach()

## Remove two outliers : Harris and Dallas
txDataR <- subset(txData, county != "Harris" & county != "Dallas")
txPanelR <- pdata.frame(txDataR, index = c("county", "year"), drop.index = TRUE, row.names = TRUE)

job.retwfe <- plm(fm,txPanelR, model = "within", effect = "twoways")
summary(job.retwfe);

job.retwre <- plm(fm,txPanelR,model = "random",random.method = "amemiya", effect = "twoways")
summary(job.retwre)

## For counties which have well drilling activities, analyse the impact of gas, oil and other wells.

## windCounty <- subset(txData, wells != 0)
## txShaleWell <- pdata.frame(txShaleWell, index = c("county", "year"), drop.index = TRUE, row.names = TRUE)
## fm.welltype <- emp ~ gas + oil

## well.twfe <- plm(fm.welltype,txShaleWell, model = "within", effect = "twoways")
## summary(well.twfe);

## well.re <- plm(fm.welltype,txShaleWell, model = "random")
## summary(well.re);

windmat = as.matrix(txPanel$windInd)
txWindCty = txData
for (i in 1:254) {
  if (all(windmat[i,] == 0)) {
   cty =  dimnames(windmat) [[1]][i]
  txWindCty = subset(txWindCty, txWindCty$county != cty)
 }
}

txWindCty <- pdata.frame(txWindCty, index = c("county", "year"), drop.index = TRUE, row.names = TRUE)

cty.twfe <- plm(fm,txWindCty, model = "within", effect = "twoways")
summary(cty.twfe);

cty.twre <- plm(fm,txWindCty, model = "random", effect = "twoways")
summary(cty.twre);

## Tests
## Tests for individual and time effects
plmtest(job.pool,effect = "twoways", type = "bp")
plmtest(job.pool,effect = "individual", type = "bp")
plmtest(job.pool,effect = "time", type = "bp")
pFtest(job.fe,job.pool)
pFtest(job.twfe,job.pool)

## Hausman test
phtest(job.fe,job.re)
phtest(job.twfe,job.twre)
phtest(cty.twfe,cty.twre)

## Tests of serial correlation
## joint/conditional LM test, test for ramdom effects and serial correlation under normality and homoskedasticity of the idiosyncratic errors
pbsytest(fm, txPanel, test = "j")
pbsytest(fm, txPanel, test = "re")
pbsytest(fm, txPanel)
## General serial correlation tests
pbltest(fm, txData, alternative = "onesided")
pbgtest(job.twfe, order = 2)
## For "short" panels with small T and large n
pwartest(fm, txData)
pwfdtest(fm, txData)
## Test for cross-sectional dependence
pcdtest(fm, txData, model = "within", effect = "twoways")
## Robust cov matrix estimation
coeftest(job.twre,vcovHC)


#######################
## ML Implementation ##
#######################

## spacial weights matrix
wmatl <- mat2listw(wmat);

## Spatial error model, random effects
## Results for the Kapoor et al. (2007)

errorre <- spml(formula = fm, data = txData, index = NULL, listw = wmatl,  model = "random", lag = F, spatial.error = "kkp")
summary(errorre)

lagre <- spml(formula = fm, data = txData, index = NULL, listw = wmatl,  model = "random", lag = T, spatial.error = "none")
summary(lagre)

## Spatial lag & error model, random effects
errlagre <- spml(formula = fm, data = txData, index = NULL, listw = wmatl,  model = "random", lag = TRUE, spatial.error = "kkp")
summary(errlagre)

## ## Spatial individual effects mu_i only, fixed effects
## spacefe <- spml(formula = fm, data = txData, listw = wmatl, model = "within", effect = "individual", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)
## summary(spacefe)

## ## Time effects gamma_t only, fixed effects
## timefe <- spml(formula = fm, data = txData, listw = wmatl, model = "within", effect = "time", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)
## summary(timefe)

## Two-way effects, fixed effects
lagtwfe <- spml(formula = fm, data = txData, index = NULL, listw = wmatl, model = "within", lag = T, spatial.error = "none", effect = "twoways", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)

errortwfe <- spml(formula = fm, data = txData, index = NULL, listw = wmatl, model = "within", lag = F, spatial.error = "kkp", effect = "twoways", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)

## Two-way spatial lag & error model, fixed effects
errlagtwofe <- spml(formula = fm, data = txData, index = NULL, listw = wmatl, model = "within", lag = TRUE, spatial.error = "b", effect = "twoways", method = "eigen", na.action = na.fail, quiet = TRUE, zero.policy = NULL, interval = NULL, tol.solve = 1e-10, control = list(), legacy = FALSE)

summary(errlagtwofe)

eff <- effects(errlagtwofe)
eff

## #######################
## ## GM Implementation ##
## #######################

## ## Spatial error model, random effects

## GM.error <- spgm(formula = fm, data = txData, listw = wmat, moments = "initial", model = "random", spatial.error = TRUE)

## summary(GM_error)

## ## #Spatial error & lag  model, random effects
## GM_full <- spgm(formula = fm, data = txData, listw = wmat, lag = TRUE, moments = "fullweights", model = "random", spatial.error = TRUE)

## summary(GM_full)

## #Spatial error model, fixed effects#
## GM_error <- spgm(formula = fm, data = txData, lag = TRUE, listw = wmat, model = "within", spatial.error = TRUE)

## summary(GM_error)

#Testing#

## LM test, no random effects assuming no spatial correlation
## H0: sigma^2_mu = 0|lambda = 0
test1 <- bsktest(x = fm, data = txData, listw = wmatl, test = "LM1")
print(class(test1))
test1

## LM test, no spatial correlation assuming no random effects
## H0: lambda = 0|sigma^2_mu = 0
test2 <- bsktest(x = fm, data = txData, listw = wmatl, test = "LM2")
test2

## Joint LM test
test3 <- bsktest(x = fm, data = txData, listw = wmatl, test = "LMH")
test3

## Conditional LM_lambda test#
## H0: lambda = 0| sigma^2_mu >= 0 
bsktest(x = fm, data = txData, listw = wmatl, test = "CLMlambda");

## Conditional LM_mu test#
## H0: sigma^2_mu = 0| lambda may or may not be zero 
bsktest(x = fm, data = txData, listw = wmatl, test = "CLMmu");


## #Spatial Hausman test# 
## testH <- sphtest(x = fm, data = txData, listw = wmatl, spatial.model = "error", method = "GM")
## testH

## mod1 <- spgm(formula = fm, data = txData, listw = wmat, lag = TRUE, moments = "fullweights", model = "random", spatial.error = TRUE)

## mod2 <- spgm(formula = fm, data = txData, listw = wmat, lag = TRUE, moments = "fullweights", model = "within", spatial.error = TRUE)

## test2 <- sphtest(x = mod1, x2 = mod2)
## test2

## #Linear hypothesis testing# 
## coeftest(errlagtwofe)

sink()
